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Abstract 

A simple 3D-reconstruction method for gamma-ray induced air showers is presented, 
which takes full advantage of the assets of a system of Atmospheric Cherenkov Tele- 
scopes combining stereoscopy and fine-grain imaging like the High Energy Stereo- 
scopic System (H.E.S.S.). The rich information collected by the cameras allows to 
select electromagnetic showers on the basis of their rotational symmetry with respect 
to the incident direction, as well as of their relatively small lateral spread. In the 
framework of a 3D-model of the shower, its main parameters — incident direction, 
shower core position on the ground, slant depth of shower maximum, average lat- 
eral spread of Cherenkov photon origins (or "photosphere 3D- width" ) and primary 
energy — are fitted to the pixel contents of the different images. For gamma-ray 
showers, the photosphere 3D- width is found to scale with the slant depth of shower 
maximum, an effect related to the variation of the Cherenkov threshold with the 
altitude; this property allows to define a dimensionless quantity uj (the "reduced 
3D-width"), which turns out to be an efficient and robust variable to discriminate 
gamma-rays from primary hadrons. In addition, the oj distribution varies only slowly 
with the gamma-ray energy and is practically independent of the zenith angle. The 
performance of the method as applied to H.E.S.S. is presented. Depending on the 
requirements imposed to reconstructed showers, the angular resolution at zenith 
varies from 0.04° to 0.1° and the spectral resolution in the same conditions from 
15% to 20%. 
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1 Introduction 



Important progress has recently been achieved in gamma-ray astronomy above 
100 GeV due to the performance of new stereoscopic systems of Atmospheric 
Cherenkov Telescopes (ACT's) equipped with high- definition imaging cam- 
eras, such as CANGAROO-3 [1], H.E.S.S.[2] and VERITAS [3]. These setups 
bring significant improvements in fiux sensitivity, angular resolution and en- 
ergy resolution, with respect to previous experiments in ground-based gamma- 
ray astronomy The increase in sensitivity is directly related to the capability 
of rejecting hadron-induced showers on the basis of the image shapes and, at 
least for known point-like sources, of shower directions. Most presently avail- 
able results from ACT's were obtained with this last constraint, the signal 
being extracted from the distribution of the pointing angle (the so-called "a 
or 6"^ plots"). However, stereoscopic systems now provide images of extended 
sources — e.g. supernova remnants [4] — and discover unexpected sources [5]; 
in such studies, hadronic showers must be rejected without using the shower 
direction. In the present article, we refer to this context as well as to that of 
point-like sources. 

The stereoscopic observation directly provides a simple geometrical recon- 
struction method based on the Hillas parameters [6] from the different im- 
ages; this method has been applied in the HEGRA experiment [7] and is 
presently adapted to H.E.S.S.[12]. High-definition imaging provides additional 
constraints which have already proved very useful, even with a single telescope, 
as in the CAT experiment [8]; in particular, the longitudinal light profile in 
the image can be modeled, providing a likelihood parameter discriminating 
gamma-rays from hadrons as well as a simultaneous determination of the 
direction of the primary 7-ray, of the impact parameter and of the shower 
energy. This method was recently extended to H.E.S.S., both in single and 
multi-telescope modes [9]. 

The method described in this article is based on a simple SD-modeling of the 
Cherenkov light emitted by an electromagnetic air shower. The rich informa- 
tion contained in several fine-grained images of such a shower provides enough 
constraints to allow an accurate reconstruction, even by means of a simple 
model incorporating the rotational symmetry of the electromagnetic cascade 
with respect to its incident direction. This allows to select 7-ray-induced show- 
ers on the basis of only two criteria with a direct physical meaning: rotational 
symmetry and small lateral spread. In section 2, the simple assumptions used 
in the shower model are described and justified. In section 3, on the basis of 
simulated gamma-rays and of real H.E.S.S. data from the blazar PKS 2155- 



* Corresponding author: M. Lemoine-Goumard, e-mail address: 
Iemoine@poly.in2p3.fr 



2 



304, it is shown that a single variable a;, directly related to the lateral spread 
of Cherenkov photon origins, is an ideal parameter for discriminating primary 
7-rays from hadrons. In section 4, the performance of the method in terms 
of gamma-ray reconstruction efficiency, of angular resolution and of hadronic 
rejection is evaluated. Finally, the spectral resolution is discussed. 



2 The reconstruction method 



2.1 Modeling assumptions 



In order to predict the distribution of Cherenkov light in the cameras of the 
different telescopes as expected from a 7-ray shower, one essentially needs 
the spatial distribution of the emission points of Cherenkov photons, and the 
angular distribution of these photons with respect to the shower axis. The most 
important characteristics of electromagnetic showers is that the distributions 
of secondary particles are on average rotationally symmetric with respect to 
the shower axis. The 3D-model of 7-ray showers presented here is based on 
the simplifying assumptions indicated below: 

(1) The emission points of Cherenkov photons are distributed according to a 
3-dimensional Gaussian law with rotational symmetry with respect to the 
shower axis, thus characterized by the following parameters (figure 1): 

• the polar angles ^0 and 0o of the shower axis in the reference frame of 
the stereoscopic system; 

• the coordinates xq and yo (in the same frame) of its intersection I with 
the ground ( "shower core position on the ground" ) ; 

• the position of the barycentre B ( "shower maximum" ) on the axis, given 
by the distance h — IB; 

• the longitudinal (o"/,) and two-fold degenerate transverse {(Tt) standard 
deviations of the Gaussian distribution, referred to as "3D-length" and 
"3D-width" respectively; 

• the total number Nc of Cherenkov photons emitted by the shower. 
The preceding quantities, referred to as "shower parameters", will be 
determined by the maximum likelihood fit described in section 2.3. 

(2) The angular distribution of Cherenkov photons, with respect to the shower 
axis, is assumed to be independent of the position of the emission point 
and of the energy of the primary 7-ray; its form, which will be given later, 
is characterized by a single parameter depending on the zenith angle only. 
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Fig. 1. electromagnetic shower in the system reference frame. 
2.2 Discussion of the preceding assumptions 

2.2.1 Assumption (1) 

The average longitudinal development of an electromagnetic shower, i.e. the 
number of electrons and positrons as a function of altitude is given by the 
Greisen formula [11] and the density profile of the atmosphere. The number of 
Cherenkov photons emitted per unit length along the shower axis follows a sim- 
ilar distribution; the variation of the Cherenkov threshold with altitude results 
in a slight shift of about 0.3 radiation lengths downwards of the Cherenkov 
profile with respect to the profile. The average longitudinal profile of a 
200 GeV vertical 7-ray shower is shown in figure 2. This profile is well de- 
scribed by a Gaussian distribution (dotted line in figure 2) for altitudes lower 
than 12000 m, particularly in the vicinity of the shower maximum which gives 
the dominant contribution to the number of collected Cherenkov photons. As 
a matter of fact, the part of the shower above 12000 m only contributes a few 
percent of collected photons. Assumption (1) is thus well justified as far as the 
longitudinal profile is concerned; the average 3D-length obtained from profiles 
such as the one shown in figure 2 is of the order of 3000 m at zenith, almost 
independent of the primary energy. 

On the other hand, the photosphere 3D-width ctt represents the average lateral 
spread of the origins of Cherenkov photons. The lateral profile is assumed to be 
Gaussian; this is not exact since the real distribution is more sharply peaked 
close to the axis; however, it will be shown (from simulations as well as from 
genuine H.E.S.S. data) that the typical 3D-width of electromagnetic showers 
at zenith is of the order of 10 m, and in most cases smaller than 15 m. Such 
a structure at a typical distance of 10 km from the telescope is viewed under 
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Fig. 2. Average longitudinal profile of 
a 200 GeV vertical 'j-ray shower: the 
logarithm of the number of electrons 
and positrons is plotted as a function of 
altitude (solid line). The approximation 
of the distribution of Ng by a Gaussian 
is shown by the dotted line. 
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Fig. 3. Probability of emission per unit 
solid angle of Cherenkov photons as 
a function of the angle e with re- 
spect to the shower axis. The distri- 
bution is shown for a shower in the 
500 GeV-1 TeV energy range at three 
different levels of its development, cor- 
responding to shower ages of 0.75, 1.0 
and 1.25. The abscissa is the ratio 
X = e/rj, where r] is the maximal 
Cherenkov angle at the altitude consid- 
ered. The graph of the function I{e) 
used in the reconstruction procedure is 
shown by the bold line. 



an angle of 1.5 mrad, smaller than the pixel size. Thus, there is no need for 
a more accurate description of the lateral distribution. The photosphere 3D- 
width, obtained with the constraint of rotational symmetry, is particularly 
useful to distinguish 7-ray from hadron showers which, in most cases, are 
much broader. 



2.2.2 Assumption (2) 

A priori, assumption (2) looks more difficult to justify, since the angular distri- 
bution of Cherenkov photons with respect to the shower axis results from that 
of electrons (which depends upon the position within the shower) combined 
with the Cherenkov cone of each electron whose half-angle depends on the 
local atmospheric density. It is thus appropriate to use the variable x — ejr]., 
in which e is the angle between the direction of the Cherenkov photon and the 
shower axis, and ri is the maximal Cherenkov angle at the altitude of the emis- 
sion point (cos?7 = l/n(2;), n{z) being the refraction index at altitude z). The 
lateral position in the shower should not be essential for the reasons explained 
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above; on the other hand, the longitudinal position can be characterized by 
the shower age, as defined in reference [11]. For a given age, one is interested 
in the probability I{e) of emission of a Cherenkov photon per unit solid angle 
around a direction making an angle e with the shower axis. The normalization 
is defined by: 

J I{e) dn^2T: j I{e) ede^l. 



Simulations show that, expressed as a function oi x — e/r), the functions / 
at fixed age do not depend strongly on the initial energy. Some of them are 
shown in figure 3 for showers in the 500 GeV-1 TeV energy range at different 
ages. The distributions are similar except for the peak at a; = 1, which rapidly 
decreases when the shower develops. Taking this peak into account in the 3D- 
model did not improve the results significantly. Therefore, we use the following 
universal function I{x): 

I{e) = K if e<7] 



K — exp 



e — f] 



if £ > 77 



in which K is fixed by the normalization: K = 1/(9711]'^). The graph of this 
function (the thick line in figure 3) is in reasonable agreement with the curves 
obtained from simulations, except in the vicinity oi e = r]. This simplified 
form actually neglects the contribution of the narrow peak around e — rj 
which is only important for small ages, i.e. for high altitudes. Ideally, rj should 
be calculated at the altitude of shower maximum, which is a parameter of the 
likelihood fit. However, coupling the angular distribution of Cherenkov light 
to the altitude often prevents the maximization procedure from converging 
and it turned out to be sufficient to use a value of 77 independent of the 
altitude and of the energy. Its dependence on the zenith angle ( was chosen 
as 77 = 15 mrady^cos^, a relation which follows if the atmospheric density is 
assumed to decrease exponentially with altitude. 



2.3 Implementation of the method 



The preceding 3D-model enables us to work out the expected number of 
Cherenkov photons qth collected by a given pixel of a given telescope as a 
function of the shower parameters listed in section 2.1. The calculation is 
done in the reference frame of the telescope of interest as shown in figure 4. 
The average direction of observation of the pixel of interest makes an angle 
9 with the telescope axis and an angle e with the shower axis. We have to 
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Fig. 4. Calculation of the expected number of collected photons in a pixel, as a 
function of the shower parameters. 

sum up the contributions of all photons emitted within the solid angle Aujpix 
covered by the pixel, provided they are pointing towards the telescope mirror 
of area Stei- At the distance r from the telescope along the line of sight, in a 
volume r'^drAcUpix, the density of Cherenkov photons ndr) is obtained from 
the shower parameters, using assumption (1). From an emission point E in 
this volume, the mirror is viewed under the solid angle dQ = Stei cos^/r^ and 
the fraction of those photons reaching the mirror is given by I{e)Stei cos^/r^. 
Integrating along the line of sight, one finds the expected value of the number 
of photons collected by the pixel of interest: 



CX) 

qth = j nc{r) r'^dr Aujpi^ I{e) 



Stei COS 9 



oo 

= Stei Aupix I{e) COS 9 j nc{r) dr 



(1) 



The last integral is easily calculated, due to assumption (1), as shown in 
Appendix 1. 



The quantities qth for all pixels, calculated for a set of shower parameters, 
are further used to build up a likelihood function for each event including at 
least two images of a given shower. Images have been previously submitted 
to a cleaning procedure in order to remove isolated clusters of pixels with a 
small charge which likely result from the night sky background. In the case 
of H.E.S.S. data analyzed in this article, the cleaning criteria were similar to 
those of reference [12] but slightly looser: pixels in the image were required to 
have a charge content above 5 photoelectrons and to have a neighbour above 
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7 photoelectrons; conversely, pixels above 7 photoelectrons were required to 
have a neighbour above 5 photoelectrons; otherwise, their charge contents 
were cleared. Pixels retained by this cleaning procedure, as well as their im- 
mediate cleared neighbours, were further used in the likelihood calculation; 
on the other hand, pixels invalidated by the calibration procedure were not 
included. The fluctuations considered in the fit are Poisson fluctuations on the 
predicted number of photons falling onto a pixel and those due to the photo- 
tube responses, the latter being considered as Gaussian. Fluctuations in the 
shower development are not taken into account, nor are correlations between 
the contents of different pixels. This has little effect on the following analysis 
since no goodness-of-fit parameter is used to discriminate gamma-ray showers 
from hadronic ones. The likelihood function for each event, described in Ap- 
pendix 2, is then maximized with respect to the 8 shower parameters defined 
in section 2.1. This is achieved by means of the MINUIT program from CERN 
[10]. The simple assumptions of the 3D-model allow a rather fast processing 
of the events. 

It should be emphasized that the present SD-reconstruction of showers is rela- 
tively insensitive to the presence of invalidated pixels in cameras; for example, 
an extreme situation in which, in all cameras, 20% of the pixels would be dam- 
aged would mainly result in a degrading of the angular resolution (an increase 
from 0.12° to 0.15° of the angular error) with only a 10% loss in reconstruction 
efficiency. 



3 Gamma-ray/ hadron discrimination based on shower shape 

3. 1 Longitudinal development: physical condition for gamma-rays 

Simulations of gamma-ray-induced showers at different energies and zenith 
angles were used to produce images in the H.E.S.S. telescopes under nor- 
mal experimental conditions (including light pollution by the night sky back- 
ground) . Those images were submitted to the cleaning procedure and analyzed 
according to the 3D-model. About 90% of the simulated showers were success- 
fully reconstructed. Since only one length scale, namely the radiation length, 
governs the development of electromagnetic showers both longitudinally and 
laterally, it is convenient to express characteristic lengths in units of radia- 
tion lengths (or equivalently in "air thickness" in g cm~^). The longitudinal 
development is thus characterized by the slant depth Dg of shower maximum 
(i.e. the thickness of air between the top of the atmosphere and the maxi- 
mum, as measured along the shower axis); it is calculated from the altitude 
Zmax of the barycentre (shower maximum) obtained by the likelihood fit and 
from the zenith angle C,. The relevance of the altitude reconstruction is illus- 
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Fig. 5. Gamma-ray showers simulated at zenith with a differential energy spectrum 
proportional to Eq'^''^ : depth of shower maximum Dg reconstructed by the fit, as a 
function of In Eq . The line shows the relation between the average value of Dg and 
IuEq given by the Greisen formula [11]. 

trated in figure 5 in which Dg is plotted as a function of the logarithm of 
the true energy £"0 for a sample of gamma-rays generated with a differential 
spectrum proportional to Eq'^''^. As expected from the Greisen formula [11], 
the average depth of shower maximum increases linearly with Ini^o- A similar 
correlation is obtained between Dg and the logarithm of the total number of 
Cherenkov photons Nc, both quantities being reconstructed by the fit; this is 
shown in figure 6 for several values of the primary gamma-ray energy Eq. This 
property allows to define an additional constraint characterizing gamma-ray 



showers: the fitted parameters Dg 
the following condition : 



(in g cm ) and Nc are required to satisfy 



61.5 (In AT^- 13.5 + T(C)) < Dg < 61.5 (In A^^ - 10 + T(C)) + 300 



(2) 



with T{() = 3.28(1 — cosC), this last term being introduced since the relation 
between the reconstructed value of Ac and the estimated value of the primary 
energy depends on (. Events satisfying condition (2) arc represented in the 
region between the two straight lines shown in figure 6 for showers at zenith. 
The fraction of gamma-rays reconstructed by the fit (about 90%) is almost 
unaffected when requiring condition (2) which results in an additional loss of 
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Fig. 6. Simulated gamma-ray showers: 
depth of shower maximum Dg as a 
function of the logarithm of the total 
number Nc of Cherenkov photons (both 
reconstructed by the fit), for different 
values of the primary energy Eq (from 
left to right: 200 GeV, 500 GeV, 1 TeV, 
5 TeV) at zenith. Condition (2) cor- 
responds to the region between the two 
straight lines; the upper bound is more 
restrictive for j-rays, due to the accu- 
mulation of background events at large 
depths. 



Fig. 7. Showers from a field of view 
free of j-ray source (far off the Galac- 
tic Plane) reconstructed as if they 
were electromagnetic showers: depth of 
shower maximum Dg as a function of 
the logarithm of the total number Nc of 
Gherenkov photons (both reconstructed 
by the fit). The observation was per- 
formed at a zenith angle of Iff'. Con- 
dition (2) corresponds to the region be- 
tween the two straight lines. 



2% for zenith angles lower than 50°. On the other hand, the fit, as applied 
to hadronic showers, sometimes converges towards very small or, more often, 
towards very large depths^ (figure 7). By requiring convergence with fitted 
parameters satisfying condition (2), one removes 70% of hadronic showers. 
The remaining ones are compatible with rotational symmetry and condition 
(2); however, their 3D- width distribution, much broader than that of gamma- 
rays, offers additional rejection criteria which we investigate now. 



3.2 Lateral development of the shower 



3.2.1 The 3D-width as a gamma-ray /hadron discriminating variable 

The potentiality of the 3D-width ctt for discriminating gamma-rays from 
hadrons is illustrated in figures 8, 9, 10, and 11, obtained from the analysis 
of H.E.S.S. data (4 hours live time) taken in 2004 on the blazar PKS2155-304 

^ This is why condition (2) is more restrictive for 7-ray showers at large depths. 
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Fig. 8. Distribution of the 3D -width from 
data taken on PKS2155-304 in 2004 . (a) 
Dotted line: source region, (h) Solid line: 
background from 5 control regions in the 
same data sample, rescaled by a factor 
1/5. 



Fig. 9. Experimental distribution of ax 
for gamma-rays (dotted line) compared 
to that obtained from simulations (solid 
line). 




100 



Fig. 10. Distribution of ax for 
gamma-rays (dotted line) obtained 
from histograms of figure 8 by subtract- 
ing the background from on-source data; 
the background distribution is given by 
the solid line. 



Fig. 11. Same histograms as in figure 10 
restricted to events with at least three 
available images. 



with a mean zenith angle of 30°. In order to monitor the hadronic background, 
data were taken with the telescope axis shifted by an angle a = 0.5° from the 
position of the source ("wobble mode"); in the following, a will be referred to 
as the offset angle. The background contribution was obtained from 5 "off" 
regions with the same acceptance as the source region. Figure 8 compares 
distributions in the source region (dotted line) and in 5 background control 
regions (solid line) rescaled by a factor 1/5. The experimental ctt distribution 
for gamma-rays, obtained by subtraction, is shown in figure 9 and found to 
be in very good agreement with that obtained from simulations based on the 
spectrum of PKS2 155-304 as measured in reference [12]. The 3D-widths of 
gamma-ray showers are essentially smaller than 20 m, thus contrasting with 
those of hadron showers, as shown in figure 10. The separation between both 
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populations is further increased by restricting the histogram to those events 
observed by at least 3 telescopes (figure 11); in this case, the corresponding 
distribution for gamma-rays is almost unaffected, whereas that of hadrons is 
shifted towards higher values. This is easily understood since gt has a real 
physical meaning for gamma-ray showers and the number of available images 
only affects the reconstruction error; on the other hand, for hadrons satisfy- 
ing the preceding fit, the assumption of rotational symmetry is generally not 
valid and the reconstructed value depends on the observation conditions. In 
any case, most hadronic showers are found to be much broader than gamma- 
ray showers. 



3.2.2 The reduced 3D-width 

Since the preceding fit also yields the altitude of shower maximum, it is con- 
venient to express ar in "air thickness" in g cm~^, i.e. to use the quantity 
a'rp = (7t p{zmax)i which p{zmax) IS the density of air at the altitude z^ax 
of the barycentre. For simulated gamma-ray showers of different primary en- 
ergies, the depth of shower maximum Dg is plotted versus in figure 12, 
both quantities being reconstructed by the fit and expressed in g cm~^. The 
average value of is found to increase with Dg. This is clearly an effect of 
the decrease of the Cherenkov threshold with depth; at larger depths, more 
low-energy electrons further from the axis produce Cherenkov light. Simula- 
tions show that, on average, a'j. scales with Dg (figure 12). Consequently, the 
dimensionless ratio u) = a'j'/Dg follows a distribution which, in the absence 
of reconstruction errors, would be independent of the gamma-ray energy and 
zenith angle. In the following, the quantity u will be referred to as the "re- 
duced 3D-width" . Figure 13 shows the ax distributions for 1 TeV gamma-ray 
showers at zenith angles 0° , 46° and 60°; in contrast, figure 14 shows the 
corresponding distributions for the reduced 3D-width a;, almost identical as 
can be seen from table 1; only a small extension towards higher values appears 
at large zenith angles. Furthermore, these distributions depend only shghtly 
on the gamma-ray energy, as shown in table 2, the average value increasing 
linearly with energy. Consequently, for a large range of primary energies, the 
uj distribution for gamma-rays is essentially concentrated in the region 

0.8 X 10"^ < < 2 X 10-2 (3) 



Of course, the upper limit is the most relevant for eliminating hadrons, the 
lower limit being nevertheless useful for events viewed by 2 telescopes only 
(see figure 28). 
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Fig. 12. Simulated gamma-ray showers: depth of shower maximum Dg as a function 
of the 3D-width a'j, (both reconstructed by the fit and expressed in g cm~'^), for 
different values of the primary energy Eq: from top to bottom, 200 GeV, 500 GeV, 
1 TeV, the last plot is the superposition of the preceding ones. The straight lines are 
intended to guide the eye. 
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Fig. 13. 3D width distributions for 
1 TeV gamma- ray showers at zenith an- 
gles 0° , 46° and 60° (from left to right). 

Table 1 

Average values and standard devia- 
tions of the CO distributions for 1 TeV 
gamma-rays. 




0.OD05 0.001 0.0015 0.002 0.0025 0.003 0.0035 0.004 0.0045 0.005 

reduced 3D-width 



Fig. 14. Distributions of the reduced 
3D-width at zenith angles 0° , 46° and 
60° for 1 TeV ^ -ray showers. 

Table 2 

Average values and standard deviations 
of the uj distributions for gamma-rays 
at zenith. 



Zenith 


< u; > 


a{u;) 


Energy 


<UJ > 


a{uj) 


angle 


xlO^ 


xlO^ 


(GeV) 


xlO^ 


xlO^ 


0° 


1.55 


0.206 


200 


1.516 


0.301 


46° 


1.57 


0.205 


500 


1.529 


0.254 


60° 


1.54 


0.212 


1000 


1.550 


0.206 



Typical statistical errors on < oo > are 
of the order of 0.01 x 10~^; those on 
a{u}) are of the order 0/O.O6 x 10~^. 



5000 



1.650 0.188 



Typical statistical errors on < lu > are 
of the order of 0.004 x 10"'^; those on 
(j{uj) are of the order of 0.03 x 10~^. 



4 Performance of the method 



4.1 The importance of the number of stereoscopic views 



The performance of the method in terms of angular or spectral resolution 
depends on the number rix of telescopes triggered by a shower; this number, 
referred to as "telescope multiplicity" in the following, depends in turn on 
the configuration of the stereoscopic system and on the position of the shower 
core on the ground. For the H.E.S.S. experiment, this is illustrated in figure 
15 showing the distributions on the ground of the impacts of vertical 7-ray 
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Fig. 15. Fraction of vertical ^ -ray showers triggering ht telescopes as a function of 
the impact position (X and Y coordinates) on the ground. From left to right: ut = 4, 
ut = 3, riT = 2. The ^-ray spectrum is a power law with a spectral index of 2.4- 
The black dots indicate the position of the H.E.S.S. telescopes. 

showers with riT = 2, = 3 and ut — ^ respectively; in this example, 7-rays 
were simulated with a power-law spectrum of differential spectral index 2.4. 
Events triggering 4 telescopes are concentrated in the central region of the 
array, whereas events triggering 2 telescopes only are peripheral and extend 
further away from the centre at higher energies. Of course, such peripheral 
events are not so accurately reconstructed as the central ones; therefore, we 
shall investigate how the angular and spectral resolution are modified when 
a minimal value of tit is required. It should be noted that the improvement 
obtained with n-r = 4 is not only due to the redundancy of the reconstruction, 
but also to the relative position of the shower core with respect to the array. 

4-2 Gamma-ray selection efficiency (shower shape) 

The similarity between the distributions of u at different zenith angles and 
energies allows to select gamma-rays by using the same criterium (3) on uj 
for all conditions of observation. On the basis of gamma-ray simulations at 
different zenith angles and energies, the selection efficiency Es for gamma-rays 
is defined as the fraction of events accepted by the fit and satisfying conditions 
(2) and (3). At this level, the selection is only based on the shower shape, not 
on its direction, as is relevant in the study of extended sources ^ . Figure 16 
shows the variation of Eg at zenith as a function of the primary energy for 
showers triggering 2, 3 and 4 telescopes. Similarly, the variation of E g as a 
function of the zenith angle is shown in figure 17 for 1 TeV 7-rays observed on 
axis; finally, for 1 TeV 7-rays observed off axis at zenith, figure 18 shows the 
variation of Sg with the offset angle, almost neghgible for showers triggering 3 

^ The global selection efficiency including the directional information (relevant for 
a point-like source) is discussed later in section 4.4. 
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Fig. 16. Reconstruction efficiency for gamma-rays at zenith, as a function of the 
primary energy, for showers triggering 2 telescopes (triangles) , 3 telescopes (squares) 
and 4 telescopes (circles) respectively. 

or 4 telescopes. The relatively high value of Eg (> 80% between 100 GeV and 
3 TeV at zenith) and its smooth variation over a large range of energies and 
zenith angles are particularly well suited for spectral analysis. 

4-3 Angular resolution 

The angular resolution is characterized by the "point-spread function" (PSF), 
i.e. the distribution of the angle 9 between the reconstructed direction and that 
of the source. It is obtained from data taken on a distant extragalactic source, 
e.g. PKS2 155-304. For the pupose of background monitoring, the source was 
observed at an angular distance a = 0.5° from the telescope axis (offset angle). 
Figure 19 shows the 6"^ distribution from the H.E.S.S. data taken on this source 
and already referred to in section 3.2.1. The angular resolution is further 
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Fig. 17. Reconstruction efficiency for 
1 TeV gamma-rays on axis as a func- 
tion of the zenith angle, for show- 
ers triggering 2 telescopes (triangles), 
3 telescopes (squares) and 4 telescopes 
(circles) respectively. 
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Fig. 18. Reconstruction efficiency for 
1 TeV gamma-rays at zenith as a 
function of the offset angle, for show- 
ers triggering 2 telescopes (triangles), 
3 telescopes (squares) and 4 telescopes 
(circles) respectively. 



improved by restricting to showers triggering at least 3 telescopes (figure 20). 
The PSF obtained from simulations in a configuration close to the one of 
PKS2155-304 (7-ray power-law spectrum with a differential spectral index of 
3.2 at 26° zenith angle) was superimposed on the histograms, showing the good 
agreement between the 6"^ distributions obtained with real data and those from 
simulations. Details of the 6"^ distribution are more clearly visible in figure 21 
in which the vertical scale is logarithmic; this distribution is obtained from 
1 TeV gamma-ray showers simulated at zenith and further reconstructed as 
explained above. It is well fitted by a linear superposition of two exponential 
laws in 6"^, 

i.e. of two Gaussian laws in projected angles {6x or 6y); the narrower one, 
referred to as "the central spot" and characterized by ai is thus superimposed 
to a broader halo characterized by (J2. The halo is more important for those 
events which trigger 2 telescopes only (figure 21 (a)); this effect, which is 
independent of the reconstruction method, is mainly related to the position 
of the impact of the shower axis on the ground, relatively far away from the 
centre of the array for those events (figure 15). On the other hand, the angular 
distribution is practically reduced to the central spot for events triggering 4 
telescopes (figure 21 (c)). Figure 22, also obtained from simulations, shows 
the variation of the spread of the central spot cti as a function of the 7-ray 
energy for different zenith angles and on-axis showers; events with all telescope 
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Fig. 19. Q'^ distribution from data 
taken on PKS21 55-304 in 2004: (a) 
Filled circles: "on source" region, (h) 
Crosses: background control regions, (c) 
The curve is the point-spread function 
obtained from simulated j-rays with a 
power-law spectrum of differential spec- 
tral index 3.2 at 26^. 



Fig. 20. 9^ distribution from data taken 
on PKS21 55-304 in 2004, restricted to 
showers triggering at least 3 telescopes: 

(a) Filled circles: "on source" region. 

(b) Crosses: background control regions. 
The curve is the point-spread function 
obtained from simulated j-rays with a 
power-law spectrum of differential spec- 
tral index 3.2 at 2(f , keeping only 
events triggering at least 3 telescopes. 



multiplicities are included. This spread is practically always smaller than 0.06° 
(about 4') and remains fairly constant at energies greater than 1 TeV. In order 
to take the halo into account, one can also characterize the angular distribution 
by the angular radius Rqs of the cone centered on the true 7-ray direction and 
containing 68% of the reconstructed axes. Its variation as a function of the 7- 
ray energy for different zenith angles is shown in figure 23 in which events with 
all telescope multiplicities are included. In the case of off-axis observations, 
the influence of the offset angle a on the angular resolution is shown in figure 
24 for 0"! and in figure 25 for Rq^, it is practically negligible for a < 1°. By 
comparing figures 22 and 23, it appears that the relative importance of the 
halo strongly depends on the zenith angle and on the position of the impact 
of the shower axis; in particular, the increase of Rqs with energy above a few 
TeV is due to the contribution of energetic showers triggering the array from 
rather remote impacts (~ 300 m from the centre). This effect is attenuated if 
one requires at least two images whose centre of gravity is within 2° from the 
camera centre, a cut used in the traditional analysis based on Hillas parameters 
(dashed hues in figure 23). 

The influence on the angular resolution of the position of the shower impact on 
the ground is shown in figures 26 and 27, obtained from simulations of vertical 
1 TeV 7-ray showers observed on-axis. Figure 26 refers to well-reconstructed 
showers, i.e. to those events whose reconstructed direction differs from the 
true one by less than 0.1°; the figure shows the fraction of this population 
among the detected events as a function of the impact position. This fraction 
is remarkably uniform — around 80% — in a large 250 m x 250 m square 
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Fig. 21. 6'^ distribution from 1 TeV j-ray showers simulated at zenith, (a) Dashed 
line: events triggering two telescopes, (b) Dashed-dotted line: events triggering three 
telescopes, (c) Solid line: events triggering four telescopes. 

area, and is still quite significant (about 40%) 300 m away from the centre 
of the array. The extension of well-reconstructed events is thus significantly 
larger than that of showers with = 4 shown in figure 27 for comparison. 



4-4 Global '-y-ray selection efficiency and hadronic rejection 



The selection criteria used in section 4.2 were based exclusively on the shower 
shape. The corresponding efficiencies averaged over a 7-ray power-law spec- 
trum with a differential spectral index of 2.4, are shown in table 3 for different 
zenith angles and for different requirements on the minimal telescope multi- 
plicity rirp. The same criteria reduce the number of hadrons by a factor Rg 
(rejection factor); the values of Rs shown in table 3 were obtained from real 
data taken in fields of view free of 7-ray sources. Whereas Ss is rather insen- 
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Fig. 22. Central spot spread ai, as a 

function of ^-ray energy for zenith an- 
gles CP (circles), ^6° (squares) and 6CP 
(triangles). 



Fig. 23. Radius of the 68% containment 

cone, as a function of ^ -ray energy for 
zenith angles (F (circles), ^6° (squares) 
and 6(P (triangles). Dotted lines show 
the effect of an additional cut requir- 
ing at least two images whose centre 
of gravity is within 2^ from the camera 
centre. 
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Fig. 24. Central spot spread ai, as 
a function of j-ray energy for offset 
angles (P (circles), 1° (squares) and 
^ (triangles); ^y-rays are simulated at 
zenith. 



Fig. 25. Radius of the 68% confidence 
level region,as a function of 'j -ray en- 
ergy for offset angles CP (circles), 1° 
(squares) and SP (triangles); 'j-rays are 
simulated at zenith. Dotted lines show 
the effect of an additional cut requir- 
ing at least two images whose centre 
of gravity is within SP from the camera 
centre. 
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Fig. 26. Fraction F of 1 TeV vertical 
7-roy showers with ^ < 0.1° as a func- 
tion of the impact position (X and Y 
coordinates in m) on the ground. 



Fig. 27. Fraction F of 1 TeV vertical 
7-ra?/ showers triggering 4 telescopes as 
a function of the impact position (X 
and Y coordinates inm) on the ground. 



sitive to TT-T, a significant improvement of tlie liadronic rejection factor Rs is 
obtained when requiring a minimal telescope multiplicity of 3. This is essen- 
tially due to the cut in the reduced 3D-width a;, since the distribution of this 
variable for hadrons strongly depends on ut as shown in figure 28; clearly, 
the constraint of rotational symmetry is more accurately checked for those 
showers surrounded by 3 or 4 telescopes. On the other hand, neither Eg, nor 
Rs vary much with the zenith angle C,. 



Table 3 

Gamma-ray selection efficiencies Eg (shape criteria) and £g (shape + direction) 
and hadronic rejection factors corresponding to the same criteria, Rg and Rg re- 
spectively, for different zenith angles and for different requirements on the minimal 
telescope multiplicity ny. Gamma-rays were simulated with a power-law spectrum 
with a differential spectral index of 2.4- Hadron rejection factors were obtained from 
real data taken in fields of view free of 'j -ray sources. Global quality factors Qg (see 
text and footnote for the definition) are also shown in the table. 



Zenith Angle 


18° 


33° 


51° 


riT 


> 2 > 3 =4 


> 2 > 3 =4 


> 2 > 3 =4 




0.86 0.90 0.90 
0.51 0.67 0.77 


0.89 0.92 0.93 
0.55 0.70 0.80 


0.85 0.89 0.90 
0.53 0.63 0.70 


Rs 

Rg 


8.9 16.1 22.4 
5000 6870 7670 


8.7 13.9 18.3 
4070 4980 4750 


8.3 11.3 13.9 
2010 2170 2880 


Qg 


36 56 67 


35 49 55 


24 30 38 
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Fig. 28. Distributions of the reduced 3D-width lo for -rays (peaks) and for hadrons 
for different observation conditions, (a) Dashed line: events triggering two tele- 
scopes, (b) Dashed-dotted line: events triggering three telescopes, (c) Solid line: 
events triggering four telescopes. Gamma-rays were simulated at 1 TeV and a zenith 
angle of 18^ ; hadrons were taken from a field of view free of -ray sources at a zenith 
angle of l(f . All histograms are normalized to the same number of events. 

In the study of a point-like source, the 7-ray selection criteria based on shower 
shape must be complemented by a cut on the angle 9 between the recon- 
structed direction and that of the source. For the sake of the present example, 
we require 9 to be lower than 0.1°/ cos^; this choice reflects the variation of 
the angular resolution with the zenith angle C,. The global 7-ray selection ef- 
ficiencies Eg and the global hadronic rejection factor i?^, obtained after this 
last cut, are shown in table 3 for different observation conditions. The global 
selection efficiency Eg is almost independent of the zenith angle and remains 
of the order of 50% with no restriction on rtx and reaches 75% for tit = 4; 
as a matter of fact, the angular cut removes most of the events from the halo 
(defined according to section 4.3) whose contribution is significant when no 
restriction on ht is applied and small for ut = 4. The global hadronic rejection 
factor varies between 2000 and 5000 with no restriction on ut (and between 
2800 and 8000 for = 4). One often defines a quality factor Qg = Eg^jR^ to 
characterize the sensitivity of an analysis method ^ . With the rather conser- 



^ The quality factor thus defined is however not a good standard of comparison 
between different experiments, since some fraction of background events is already 
removed at the trigger level. In particular, muons, which make a significant part of 
the background events for a single telescope are eliminated by the central trigger in 
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Fig. 29. Average value of In Nc as a 
function of Eq for j-ray showers simu- 
lated on axis at different zenith angles: 
0° (circles), 46° (squares) and 60° (tri- 
angles), and triggering 4 telescopes. 
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Fig. 30. Average value of In Nc as a 
function of Eq for j-ray showers sim- 
ulated at zenith for two offset angles: 
0° (circles), 1° (squares) and2° (trian- 
gles), and triggering 4 telescopes. 



vative cuts used above, Qg is of the order of 20 to 40 with no restriction on 
riT and is practically doubled when one requires rix = 4. It should be noted 
that a minimal set of cuts have been used in the preceding example and that 
no optimization procedure has been applied. Additional requirements (e.g. on 
the minimum number of photo-electrons per telescope) and different sets of 
cuts can be used to improve the analysis, depending on the intensity of the 
source. 



4-5 Gamma-ray energy measurement 



The energy Eq of the primary gamma-ray is reconstructed calorimetrically 
from the total number of Cherenkov photons obtained from the fit. If the 
3D-model of the shower described above were a perfect representation of the 
electromagnetic shower, would be, on average, almost^ proportional to 
Eq. As a matter of fact, simulations of 7-ray showers for fixed values of Eq, 
show that the average value of In Nc (noted as (InA'^c)), obtained from the 
likelihood fit, varies almost linearly with IylEq for fixed observing conditions, 
namely: telescopes multiplicity ut, direction of the telescope axes (i.e. zenith 
angle (), direction of the shower axis (i.e. offset angle a). For those showers 
triggering 4 telescopes, i.e. those for which the sampling is more homogeneous, 
the slope a = d{lnNc)/dlnEQ is actually close to 1; this is shown in figure 29 
for showers generated on-axis at different zenith angles, as well as in figure 30 
for showers at zenith generated off-axis at different angles. 



a stereoscopic system. 

The variation of the Cherenkov threshold as a function of the altitude makes this 
statement only approximate. 
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Fig. 31. Average value of In Nc as a 
function of Eq for j-ray showers simu- 
lated on axis at different zenith angles: 
0° (circles), 46° (squares) and 60° (tri- 
angles), and triggering 2 telescopes. 
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Fig. 32. Average value of In Nc as a 
function of Eq for j-ray showers sim- 
ulated at zenith for two offset angles: 
0° (circles), 1° (squares) and2° (trian- 
gles), and triggering 2 telescopes. 
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Fig. 33. Average value of In Nc as a 
function of Eq for j-ray showers simu- 
lated on axis at different zenith angles: 
0° (circles), 46° (squares) and 60° (tri- 
angles), and triggering > 3 telescopes. 
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Fig. 34. Average value of In Nc as a 
function of Eq for j-ray showers sim- 
ulated at zenith for two offset angles: 
0° (circles), 1° (squares) and 2° (trian- 
gles), and triggering > 3 telescopes. 



However, as an effect of the approximations of tfie 3D-model used in tfie fit, 
tlie coefficients a and b in tlie formula ( In A^^) = clIhEq -\- b depend on tlie 
observing conditions; tlie slight deviation from linearity at the highest energies 
for large offset angles is due to shower images not fully contained in the field 
of view. From figures 29 and 30, one verifies that the slope a is slightly higher 
than 1 for showers triggering 2 or 3 telescopes (figures 31, 32, 33, 34). For a 
given value of n^, a is almost independent of the zenith angle and of the offset 
angle. 



The impact parameter dT, i.e. the distance between the centre of the array 
and the shower axis, also affects the sampling of the shower. In figures 29 
to 34, InA'^c was averaged over all values of compatible with the trigger 
conditions; a better energy estimate is obtained if the measured value of dx is 
taken into account. For fixed values of ht, C Q^, the distributions of IniVc 
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Fig. 35. Distribution of In Nc versus the impact parameter dx for 500 GeV j-ray 
showers simulated on axis at zenith; from left to right: ut = 4:, tit = 3, tit = 2. 

are plotted for different values of dj^ in figures 35 (500 GeV on-axis showers at 
zenith); in most cases, (InA^^) varies only slowly with d^- However, at large 
distance this variation can be more important as one can see on the distri- 
bution obtained for showers triggering only 2 telescopes, as shown in figure 35. 

The energy estimator is based essentially on the value of Nc, but also on tit, (, 
a and d^. Simulations for fixed values of these parameters provide a relation 
between (InA^^c) and Inii^o which can be generalized to arbitrary values of (, 
a and dr through interpolations. This same relation is then used for a given 
event satisfying the likelihood fit to find the energy estimator Ej. as a function 
of Nc, riT, C, Oi and dr- The distribution oi\\i{Erl Eq) is found to be Gaussian 
to a good approximation for all observing conditions as shown in figure 36. 
These Gaussian distributions are characterized by their bias 5 = {\n{E^/EQ)) 
and their standard deviation a{hi{Er/ Eq)) k, AEr/Er- The variations of 6 
and of AEr/Er with energy at different zenith angles are shown in figure 37 
for on-axis showers; figure 38 shows the corresponding variations at different 
offset angles for showers at zenith (restricting to 7-ray showers with ^ < 0.1° 
and an angular distance from the centre of the camera lower than 2°). The 
bias is smaller than 5% and AEr/Er is of the order of 15 to 20% at low zenith 
angles. 



5 Conclusion 

The reconstruction method described above differs from the more traditional 
ones (ref. [12] and [9]) in several aspects: 

• The analysis is based on shower parameters in 3 dimensions, not on image 
parameters; in this way correlations between different stereoscopic views of 
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Fig. 36. Distribution ofln^Er/Eo) for 'j- ray showers simulated on-axis at zenith for 
different primary energies: 200 GeV (top left), 500 GeV (top right), 1 TeV (bottom 
left), 10 TeV (bottom right) 



the same shower are taken into account. 

• Gamma-ray candidates are selected on the basis of a few criteria based on 
physical properties: rotational symmetry of electromagnetic shower, depth 
of shower maximum, lateral spread of the Cherenkov photosphere at shower 
maximum. 

• The distribution of the ratio of the last two variables, namely the "reduced 
3D-width" LJ, is found to be almost independent of the 7-ray energy and 
zenith angle and 7-rays can be efficiently selected by requiring 0.8 x 10~^ < 
uj < 2 X 10"^ for all observing conditions. This critcrium provides a 7- 
ray/hadron discrimination based on the shower shape (thus relevant for the 
study of extended sources) and completely independent of simulations. 
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Fig. 37. Energy measurement at different zenith angles C and on-axis showers. Open 
symbols: bias 6 in ln(£'.r/£'o) a function of the true primary energy Eq; filled 
symbols: standard deviation ofln{Er/EQ) as a function of Eq. Circles are for C, = 0°, 
squares for C, = 46° and triangles for C, = 60° . 

With a minimal set of natural cuts, the selection efficiency for 7-rays as applied 
to H.E.S.S. data is rather uniform between 100 GeV and 10 TeV and of the 
order of 80% (shape criteria) and of 50% (shape and direction criteria); the an- 
gular resolution (from 0.04° to 0.1° at zenith) and the energy resolution (from 
15% to 20% in the same conditions) are comparable to those obtained with 
the standard H.E.S.S. analysis [12]. This is illustrated in a recent article on 
the blazar H2356-309 [13] where both analyses were used: using the 3D-model, 
an excess of 715 7-ray candidates is obtained with a significance of 10.9 cr, to 
be compared with 591 7-ray candidates and a significance of 9.7 a with the 
H.E.S.S. standard analysis. Furthermore, requiring at least 3 triggering tele- 
scopes, the 3D-model yields a higher significance (11.6 a) while keeping 453 
7-ray candidates. Thus, the reconstruction method explained above illustrates 
how the quality of the stereoscopy, characterized by the telescope multiphcity, 
improves the angular resolution and the hadronic rejection of the array, thus 
its sensitivity. 
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Fig. 38. Energy measurement at different offset angles a and showers at zenith. 
Open symbols: bias 6 in \n(Er/EQ) as a function of the true primary energy Eq; 
filled symbols: standard deviation ofln{Er/EQ) as a function of Eq. Circles are for 
a = 0°, squares for a = 1° and triangles for a = 2° . 
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6 Appendix 1: Calculation of the pixel content expected from the 
3D-model 



The expected number of Cherenkov photons cohected by a given pixel is given 
by formula (1) in which the notations are defined in figure 4. The formula 
includes the integral ndr) dr . taken along the line of sight corresponding 
to the pixel of interest. The calculation of this integral is explicited below. 
Let p be the unit vector along the line of sight and s the unit vector along 

— * 

the shower axis, both directed upwards. Let xb = OB be the vector defined 
by the optical centre O of the telescope and the barycentre B of the shower. 
The parameters s, a^, and Xb are provided by the 3D-model. Due to 
assumption (1) (see section 2.2), for a shower with Nc Cherenkov photons, the 
photon density is given by: 



(27r)3/Via2 V 2 



M = \ + and 



^ = {rp — Xb) ■ s — rcose — xb ■ s and 

r/^ = {rp — xbY — (rcose — xb • s)'^ ■ 

Defining now the following quantities, independent of r: 

Bs^xb-s ; Bp^ xb-P ; u^cose ; A| = - 5^ , 

2 22,2 /1 2\ 2 2i2-2 2 2 2 

(T„ = arpU + ai^[l — u ) — arp cos s + ai^sm e ; ajj — aj^ — aj^ , 



M, a second degree polynomial in r, can be written in the canonical form: 



M 



2^2 



aifcr, 



alBp - uBsgI, 



+ R, 



R being a quantity independent of r given by: 
1 



R 



KK-^U^Bp-B, 



Finally, equation (4) takes the form: 

7Veexp[-i?/2) 

f^c{r) = /^_xs/9_ _9, exp 



(27r)3/2cri4 



alBp - uB.^af ^ 



r — 
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and the integral in formula (1) is given by: 



/ nc{r) dr = J^"^ exp 



A2 



a 



cr- 



(5) 



in which C — 1 — freq ^— , the function freq(x) being defined 

1 ^ 

freq(a;) = ^= / exp{-t'^/2) dt . 
v27r J 

— oo 



as: 



7 Appendix 2: Calculation of the likelihood function 



Since correlations between pixel contents are not taken into account, the likeli- 
hood function for an event takes the form C = Yl^ii in which ii is the likelihood 
function of the individual pixel i. The product is taken over all the selected 
pixels of the different images of a given shower. For a given pixel for which 
the preceding 3D- model predicts an average number of photoelectron q^h, the 
probability to measure a charge (expressed in number of photoelectrons) in 
the interval [q, q + dq] is given by: 

^ exp{-qth) g^fe 1 



In this formula, Poissonian fluctuations of the effective number n of photoelec- 
trons have been assumed, as well as Gaussian fluctuations in the phototube 
response with a standard deviation a. The corresponding probability distribu- 
tion function at fixed q depends (through qth) on the shower parameters to be 
fitted and thus represents the factor £ of the likelihood function corresponding 
to the pixel of interest. In practice, since in H.E.S.S. a — 0.4 photoelec- 
trons, the sum over n is restricted to the integer values within a ±3(7 interval 
around q. 
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